SlingShot_Wnt11_NPs
Mehran Piran
6/27/2021
\(~\)
\(~\)
library(Seurat)## Attaching SeuratObjectlibrary(grid)
so = readRDS("/mnt/market/anclab-rstudio-server/home/swil0005/projects/Wnt11/raw_data/Wnt11/wnt11-QC-clustered-sctransform.RDS")
DimPlot(so , reduction = "umap" , label = T)
grid.text("(Fig1)", 0.15, 0.85)
DimPlot(so , reduction = "umap" , group.by = "genotype")
grid.text("(Fig2)", 0.15, 0.85)
head(Idents(so) , 20)## TGACTCCCATCATCTT AGGTTACCAGTCTGGC CTACTATCACATATGC CTCATGCTCTAACGCA
## 0 0 0 18
## TGAGGTTCAATAACGA GTGATGTGTAGTCCTA CACATGACATCTCAAG GTCACTCAGTTACTCG
## 21 2 4 3
## CACGTTCGTTCTCCCA TCATGCCCATCGCTCT CTGAGGCCAGCAGTGA TCTCTGGTCTTCCCGA
## 1 1 1 9
## CGACAGCGTAACCCGC CATAAGCCAGCGAACA GAGACCCAGTCGCCAC CGTAAGTGTAAGCAAT
## 0 5 0 0
## GAGTGTTAGCCTTCTC AAAGTCCGTAGATCCT CGAGGAACACAAGTTC GGTGTTAAGCGCACAA
## 13 8 19 2
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21head(so[[]])## orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
## TGACTCCCATCATCTT wnt11mouse 22788 4784 3072 5
## AGGTTACCAGTCTGGC wnt11mouse 15595 3744 2669 5
## CTACTATCACATATGC wnt11mouse 35180 6039 5276 5
## CTCATGCTCTAACGCA wnt11mouse 20594 4508 3863 5
## TGAGGTTCAATAACGA wnt11mouse 30360 5296 3926 5
## GTGATGTGTAGTCCTA wnt11mouse 11949 3254 1532 5
## HTO_maxID HTO_secondID HTO_margin
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT TS-A-55-GTTGTGAGCACGAGA 1.7065694
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC TS-A-54-TATTTCCACCCGGTC 1.6277081
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC TS-A-52-GTCCGACTAATAGCT 2.2671312
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC TS-A-55-GTTGTGAGCACGAGA 2.0316834
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC TS-A-53-CAGGTTGTTGTCATT 1.8152870
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC TS-A-55-GTTGTGAGCACGAGA 0.8745784
## HTO_classification HTO_classification.global
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT Singlet
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC Singlet
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC Singlet
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC Singlet
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC Singlet
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC Singlet
## hash.ID percent.mito S.Score G2M.Score
## TGACTCCCATCATCTT TS-A-53-CAGGTTGTTGTCATT 1.505178 -0.04425451 0.3310395
## AGGTTACCAGTCTGGC TS-A-51-AACCTTTGCCACTGC 2.378968 0.04743516 -0.2732218
## CTACTATCACATATGC TS-A-54-TATTTCCACCCGGTC 2.151791 -0.08918952 0.3582894
## CTCATGCTCTAACGCA TS-A-51-AACCTTTGCCACTGC 3.700107 -0.37091462 -0.3498180
## TGAGGTTCAATAACGA TS-A-51-AACCTTTGCCACTGC 2.292490 -0.27135004 0.5070627
## GTGATGTGTAGTCCTA TS-A-54-TATTTCCACCCGGTC 2.426981 -0.26513688 0.3523543
## Phase high.mito nCount_SCT nFeature_SCT SCT_snn_res.0.2
## TGACTCCCATCATCTT G2M FALSE 16238 4714 2
## AGGTTACCAGTCTGGC S FALSE 15480 3743 0
## CTACTATCACATATGC G2M FALSE 16372 5006 0
## CTCATGCTCTAACGCA G1 FALSE 16504 4448 7
## TGAGGTTCAATAACGA G2M FALSE 16269 4731 1
## GTGATGTGTAGTCCTA G2M FALSE 14626 3253 2
## SCT_snn_res.0.3 SCT_snn_res.0.4 SCT_snn_res.0.5
## TGACTCCCATCATCTT 1 1 1
## AGGTTACCAGTCTGGC 0 0 0
## CTACTATCACATATGC 0 0 0
## CTCATGCTCTAACGCA 8 8 8
## TGAGGTTCAATAACGA 2 2 2
## GTGATGTGTAGTCCTA 1 1 1
## SCT_snn_res.0.6 SCT_snn_res.0.7 SCT_snn_res.0.8
## TGACTCCCATCATCTT 0 1 0
## AGGTTACCAGTCTGGC 0 0 0
## CTACTATCACATATGC 0 0 0
## CTCATGCTCTAACGCA 12 11 18
## TGAGGTTCAATAACGA 1 19 21
## GTGATGTGTAGTCCTA 2 1 2
## SCT_snn_res.0.9 SCT_snn_res.1 SCT_snn_res.1.1 SCT_snn_res.1.2
## TGACTCCCATCATCTT 1 7 7 6
## AGGTTACCAGTCTGGC 0 7 7 6
## CTACTATCACATATGC 0 3 4 2
## CTCATGCTCTAACGCA 20 20 21 24
## TGAGGTTCAATAACGA 15 16 16 19
## GTGATGTGTAGTCCTA 1 0 0 10
## seurat_clusters genotype
## TGACTCCCATCATCTT 6 wild-type
## AGGTTACCAGTCTGGC 6 wild-type
## CTACTATCACATATGC 2 wnt11-hom
## CTCATGCTCTAACGCA 24 wild-type
## TGAGGTTCAATAACGA 19 wild-type
## GTGATGTGTAGTCCTA 10 wnt11-homso@assays## $RNA
## Assay data with 31053 features for 8065 cells
## Top 10 variable features:
## Hbb-y, Hba-x, S100a8, Fabp4, Stfa1, Cd74, Fxyd2, S100a9, Hmox1, Ifitm1
##
## $HTO
## Assay data with 5 features for 8065 cells
## First 5 features:
## TS-A-51-AACCTTTGCCACTGC, TS-A-52-GTCCGACTAATAGCT,
## TS-A-53-CAGGTTGTTGTCATT, TS-A-54-TATTTCCACCCGGTC,
## TS-A-55-GTTGTGAGCACGAGA
##
## $SCT
## Assay data with 18549 features for 8065 cells
## Top 10 variable features:
## Apoe, C1qb, Hbb-bs, Hba-a1, C1qa, Fcer1g, Lyz2, Hba-a2, Pf4, Tyrobpso@reductions## $pca
## A dimensional reduction object with key PC_
## Number of dimensions: 50
## Projected dimensional reduction calculated: FALSE
## Jackstraw run: FALSE
## Computed using assay: SCT
##
## $umap
## A dimensional reduction object with key UMAP_
## Number of dimensions: 2
## Projected dimensional reduction calculated: FALSE
## Jackstraw run: FALSE
## Computed using assay: SCTDefaultAssay(so)## [1] "SCT"# Number of unique groups in each meta.data column
lapply(apply(so[[]] , 2 , unique) , length)## $orig.ident
## [1] 1
##
## $nCount_RNA
## [1] 6897
##
## $nFeature_RNA
## [1] 3580
##
## $nCount_HTO
## [1] 3820
##
## $nFeature_HTO
## [1] 1
##
## $HTO_maxID
## [1] 5
##
## $HTO_secondID
## [1] 5
##
## $HTO_margin
## [1] 8060
##
## $HTO_classification
## [1] 5
##
## $HTO_classification.global
## [1] 1
##
## $hash.ID
## [1] 5
##
## $percent.mito
## [1] 8039
##
## $S.Score
## [1] 8065
##
## $G2M.Score
## [1] 8060
##
## $Phase
## [1] 3
##
## $high.mito
## [1] 1
##
## $nCount_SCT
## [1] 2142
##
## $nFeature_SCT
## [1] 2478
##
## $SCT_snn_res.0.2
## [1] 12
##
## $SCT_snn_res.0.3
## [1] 15
##
## $SCT_snn_res.0.4
## [1] 15
##
## $SCT_snn_res.0.5
## [1] 16
##
## $SCT_snn_res.0.6
## [1] 19
##
## $SCT_snn_res.0.7
## [1] 21
##
## $SCT_snn_res.0.8
## [1] 22
##
## $SCT_snn_res.0.9
## [1] 23
##
## $SCT_snn_res.1
## [1] 25
##
## $SCT_snn_res.1.1
## [1] 25
##
## $SCT_snn_res.1.2
## [1] 28
##
## $seurat_clusters
## [1] 28
##
## $genotype
## [1] 2which(unlist(lapply(apply(so[[]] , 2 , unique) , length)) == 22)## SCT_snn_res.0.8
## 25unique(so[[]][ , 25])## [1] 0 18 21 2 4 3 1 9 5 13 8 19 7 12 11 6 20 10 16 14 17 15
## Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21so.sub = subset(so , idents = c("1" , "5" , "13" , "15" , "18" , "21"))
so.sub## An object of class Seurat
## 49607 features across 1895 samples within 3 assays
## Active assay: SCT (18549 features, 3000 variable features)
## 2 other assays present: RNA, HTO
## 2 dimensional reductions calculated: pca, umapDimPlot(so.sub , reduction = "umap" , label = T)
grid.text("(Fig3)", 0.15, 0.85)
DimPlot(so.sub , reduction = "umap" , group.by = "genotype")
grid.text("(Fig4)", 0.15, 0.85)
DimPlot(so.sub , reduction = "umap" , group.by = "Phase")
grid.text("(Fig5)", 0.15, 0.85)
\(~\)
\(~\)
plot = DimPlot(so.sub , reduction = "umap" , label = T)
selected = CellSelector(plot = plot)
selected
so.small = CellSelector(plot = plot , object = so.sub , ident = 'SelectedCells')
head(colnames(so.small[[]]))
DimPlot(so.small)
unique(so.small[[]][,25])
lapply(lapply(apply(so.small[[]] , 2 , unique) , function(x) x %in% "SelectedCells") , which)
Idents(so.sub)
Idents(so.small)
levels(so.sub)
levels(droplevels(Idents(so.sub)))
levels(so.small)
levels(droplevels(Idents(so.small)))
so.small = subset(so.small , idents = "SelectedCells")
DimPlot(so.small)
Idents(so.small) = "SCT_snn_res.0.8"
DimPlot(so.small , label = T , label.size = 6)\(~\)
\(~\)
library(Seurat)
library(SingleCellExperiment)
library(slingshot)
library(RColorBrewer)
library(bioc2020trajectories)## Warning: replacing previous import 'HDF5Array::path' by 'igraph::path' when
## loading 'bioc2020trajectories'## Warning: replacing previous import 'igraph::compose' by 'purrr::compose' when
## loading 'bioc2020trajectories'## Warning: replacing previous import 'igraph::simplify' by 'purrr::simplify' when
## loading 'bioc2020trajectories'library(ggbeeswarm)
library(ggthemes)
library(BiocParallel)
library(tradeSeq)
library(future)
library(tidyverse)
so.small = readRDS("/mnt/market/anclab-rstudio-server/home/mpir0002/Slingshot/so.small.rds")\(~\)
\(~\)
Trajectories on Clusters 1, 5, 13, 15, 18, and 21
DimPlot(so.small , reduction = "umap" , label = T , label.size = 6)
grid.text("(Fig6)", 0.15, 0.85)
sce <- as.SingleCellExperiment(so.small)
sce## class: SingleCellExperiment
## dim: 18549 1847
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
## AGACAAAAGTACTGTC
## colData names(32): orig.ident nCount_RNA ... genotype ident
## reducedDimNames(2): PCA UMAP
## altExpNames(0):head(Idents(so.small))## CTCATGCTCTAACGCA TGAGGTTCAATAACGA CACGTTCGTTCTCCCA TCATGCCCATCGCTCT
## 18 21 1 1
## CTGAGGCCAGCAGTGA CATAAGCCAGCGAACA
## 1 5
## Levels: 1 5 13 15 18 21slingshot.1to18 <- slingshot(sce, reducedDim = 'UMAP', clusterLabels = sce$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")## Using full covariance matrixslingshot.1to18## class: SingleCellExperiment
## dim: 18549 1847
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
## AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
## slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):sds.1to18 <- SlingshotDataSet(slingshot.1to18)
sds.1to18## class: SlingshotDataSet
##
## Samples Dimensions
## 1847 2
##
## lineages: 2
## Lineage1: 1 21 5 13 18
## Lineage2: 1 21 5 13 15
##
## curves: 2
## Curve1: Length: 10.129 Samples: 1672.77
## Curve2: Length: 10.861 Samples: 1710.37slingshot.1to15 <- slingshot(sce, reducedDim = 'UMAP', clusterLabels = sce$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "15")## Using full covariance matrixslingshot.1to15## class: SingleCellExperiment
## dim: 18549 1847
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(1847): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... ATTCACTGTCCAGAAG
## AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
## slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):sds.1to15 <- SlingshotDataSet(slingshot.1to15)
sds.1to15## class: SlingshotDataSet
##
## Samples Dimensions
## 1847 2
##
## lineages: 2
## Lineage1: 1 21 5 13 18
## Lineage2: 1 21 5 13 15
##
## curves: 2
## Curve1: Length: 10.129 Samples: 1672.77
## Curve2: Length: 10.861 Samples: 1710.37\(~\)
\(~\)
Plotting the trajectories
colors <- rainbow(100, alpha = 1)
plot(reducedDims(slingshot.1to18)$UMAP,
col = colors[cut(slingshot.1to18$slingPseudotime_1,breaks=100)],
pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig7)", 0.15, 0.85)
plot(reducedDims(slingshot.1to18)$UMAP,
col = colors[cut(slingshot.1to18$slingPseudotime_2,breaks=100)],
pch=16, asp = 1, main = "colored_basedon_Pseudotime2_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig8)", 0.15, 0.85)
# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to18)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to18$SCT_snn_res.0.8)], pch=16,
main = "colored_basedon_Clusters_cluster1to18")
lines(SlingshotDataSet(slingshot.1to18), lwd=2)
grid.text("(Fig9)", 0.15, 0.85)
\(~\)
\(~\)
colors <- rainbow(100, alpha = 1)
plot(reducedDims(slingshot.1to15)$UMAP,
col = colors[cut(slingshot.1to15$slingPseudotime_1,breaks=100)],
pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig10)", 0.15, 0.85)
plot(reducedDims(slingshot.1to15)$UMAP,
col = colors[cut(slingshot.1to15$slingPseudotime_2,breaks=100)],
pch=16, asp = 1, main = "colored_basedon_Pseudotime2_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig11)", 0.15, 0.85)
# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to15)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to15$SCT_snn_res.0.8)], pch=16,
main = "colored_basedon_Clusters_cluster1to15")
lines(SlingshotDataSet(slingshot.1to15), lwd=2)
grid.text("(Fig12)", 0.15, 0.85)
\(~\)
\(~\)
\(~\)
\(~\)
Trajectories on Clusters 1, 5, 13, and 21
so.small2 = subset(so.small , idents = c("1" , "5" , "13" , "21"))
DimPlot(so.small2 , reduction = "umap" , label = T , label.size = 6)
grid.text("(Fig13)", 0.15, 0.85)
sce2 = as.SingleCellExperiment(so.small2)
slingshot.1to13 = slingshot(sce2 , reducedDim = 'UMAP', clusterLabels = sce2$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "13")## Using full covariance matrixslingshot.1to13## class: SingleCellExperiment
## dim: 18549 1650
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(1650): TGAGGTTCAATAACGA CACGTTCGTTCTCCCA ... ATTCACTGTCCAGAAG
## AGACAAAAGTACTGTC
## colData names(34): orig.ident nCount_RNA ... slingClusters
## slingPseudotime_1
## reducedDimNames(2): PCA UMAP
## altExpNames(0):sds.1to13 = SlingshotDataSet(slingshot.1to13)
sds.1to13## class: SlingshotDataSet
##
## Samples Dimensions
## 1650 2
##
## lineages: 1
## Lineage1: 1 21 5 13
##
## curves: 1
## Curve1: Length: 9.1458 Samples: 1650\(~\)
\(~\)
Plotting the trajectories
colors <- rainbow(100, alpha = 1)
plot(reducedDims(slingshot.1to13)$UMAP,
col = colors[cut(slingshot.1to13$slingPseudotime_1,breaks=100)],
pch=16, asp = 1, main = "colored_basedon_Pseudotime1_cluster1to13")
lines(SlingshotDataSet(slingshot.1to13), lwd=2)
grid.text("(Fig14)", 0.15, 0.75)
# Show the cells coloured by clusters
plot(reducedDims(slingshot.1to13)$UMAP, col = brewer.pal(9,'Set1')[droplevels(slingshot.1to13$SCT_snn_res.0.8)], pch=16,
main = "colored_basedon_Clusters_cluster1to13")
lines(SlingshotDataSet(slingshot.1to13), lwd=2)
grid.text("(Fig15)", 0.15, 0.75)
\(~\)
\(~\)
Compute an imbalance score to show whether nearby cells have the same condition of not
scores = imbalance_score(rd = reducedDims(sce)$UMAP, cl = colData(sce)$genotype,
k = 20, smooth = 40)
imbalance_data = merge(as.data.frame(reducedDims(sce)$UMAP), as.data.frame(scores), by = 0, all = TRUE)
ggplot() + geom_point(aes(x = UMAP_1,y = UMAP_2,colour = scaled_scores),data=imbalance_data, size = 3) + scale_colour_viridis_c(option = "plasma") + theme_linedraw() + ggtitle("Imbalance of Wild-type and Wnt11KO cells across UMAP space")
grid.text("(Fig16)", 0.15, 0.85)
\(~\)
\(~\)
suppressPackageStartupMessages({
library(scales)
library(viridis); library(UpSetR)
library(pheatmap); library(msigdbr)
library(fgsea); library(knitr)
library(ggplot2); library(gridExtra)
library(tradeSeq)
})\(~\)
\(~\)
Distribution of cells based on cell-cyle and genotype
shuffle <- sample(ncol(sce))
layout(matrix(1:2, nrow = 1))
par(mar = c(4.5,4,1,1))
plot(reducedDims(sce)$UMAP[shuffle, ],
asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2",
col = alpha(c(1:2)[factor(colData(sce)$genotype)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 1:2, bty = "n",
legend = levels(factor(colData(sce)$genotype)))
grid.text("(Fig17)", 0.15, 0.75)
plot(reducedDims(sce)$UMAP[shuffle, ], asp = 1, pch = 16, xlab = "UMAP-1", ylab = "UMAP-2",
col = alpha(c(3,4,5)[factor(colData(sce)$Phase)][shuffle], alpha = .5))
legend("topright", pch = 16, col = 3:5, bty = "n", legend = levels(factor(colData(sce)$Phase)))
grid.text("(Fig18)", 0.25, 0.75)
\(~\)
\(~\)
Plot the distribution of cells per cluster across pseudotimes.
slingshot.1to18$SCT_snn_res.0.8 = droplevels(slingshot.1to18$SCT_snn_res.0.8)
slingshot.1to15$SCT_snn_res.0.8 = droplevels(slingshot.1to15$SCT_snn_res.0.8)
slingshot.1to13$SCT_snn_res.0.8 = droplevels(slingshot.1to13$SCT_snn_res.0.8)
ggplot(as.data.frame(colData(slingshot.1to18)),
aes(x = slingshot.1to18$slingPseudotime_1,
y = forcats::fct_relevel(slingshot.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
geom_quasirandom(groupOnX = F) +
scale_color_tableau() + theme_classic() +
xlab("Slingshot pseudotime") + ylab("Clusters") +
ggtitle("Aggregated") + NoLegend()## Warning: Removed 111 rows containing missing values (position_quasirandom).grid.text("(Fig19)", 0.15, 0.85)
summary(slingshot.1to18$slingPseudotime_1)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 1.210 3.002 3.540 5.427 10.129 111so.wild = subset(so.small , subset = genotype == "wild-type")
so.wild## An object of class Seurat
## 49607 features across 827 samples within 3 assays
## Active assay: SCT (18549 features, 3000 variable features)
## 2 other assays present: RNA, HTO
## 2 dimensional reductions calculated: pca, umapso.wnt11 = subset(so.small , subset = genotype == "wnt11-hom")
so.wnt11## An object of class Seurat
## 49607 features across 1020 samples within 3 assays
## Active assay: SCT (18549 features, 3000 variable features)
## 2 other assays present: RNA, HTO
## 2 dimensional reductions calculated: pca, umapsce.wild = as.SingleCellExperiment(so.wild)
sce.wnt11 = as.SingleCellExperiment(so.wnt11)
slingshot.wild.1to18 <- slingshot(sce.wild, reducedDim = 'UMAP', clusterLabels = sce.wild$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")## Using full covariance matrixslingshot.wild.1to18## class: SingleCellExperiment
## dim: 18549 827
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(827): CTCATGCTCTAACGCA TGAGGTTCAATAACGA ... GGGCTCAAGCATCTTG
## CGAATTGGTGGATCGA
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
## slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):slingshot.wnt11.1to18 <- slingshot(sce.wnt11, reducedDim = 'UMAP', clusterLabels = sce.wnt11$SCT_snn_res.0.8 , start.clus = "1" , end.clus = "18")## Using full covariance matrixslingshot.wnt11.1to18## class: SingleCellExperiment
## dim: 18549 1020
## metadata(0):
## assays(2): counts logcounts
## rownames(18549): Xkr4 Sox17 ... CAAA01118383.1 CAAA01147332.1
## rowData names(6): sct.detection_rate sct.gmean ...
## sct.residual_variance sct.variable
## colnames(1020): CATAAGCCAGCGAACA GAGTGTTAGCCTTCTC ... ATTCACTGTCCAGAAG
## AGACAAAAGTACTGTC
## colData names(35): orig.ident nCount_RNA ... slingPseudotime_1
## slingPseudotime_2
## reducedDimNames(2): PCA UMAP
## altExpNames(0):ggplot(as.data.frame(colData(slingshot.wild.1to18)),
aes(x = slingshot.wild.1to18$slingPseudotime_1,
y = forcats::fct_relevel(slingshot.wild.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
geom_quasirandom(groupOnX = F) +
scale_color_tableau() + theme_classic() +
xlab("Pseudotime1") + ylab("Clusters") +
ggtitle("Wild Type") + NoLegend()## Warning: Removed 35 rows containing missing values (position_quasirandom).grid.text("(Fig20)", 0.15, 0.85)
ggplot(as.data.frame(colData(slingshot.wnt11.1to18)),
aes(x = slingshot.wnt11.1to18$slingPseudotime_1,
y = forcats::fct_relevel(slingshot.wnt11.1to18$SCT_snn_res.0.8,"1", "5", "13", "15", "18" , "21"), colour = SCT_snn_res.0.8)) +
geom_quasirandom(groupOnX = F) +
scale_color_tableau() + theme_classic() +
xlab("Pseudotime1") + ylab("Clusters") +
ggtitle("Wnt11KO") + NoLegend()## Warning: Removed 80 rows containing missing values (position_quasirandom).grid.text("(Fig21)", 0.15, 0.85)
ggplot(as.data.frame(colData(slingshot.1to18)),
aes(x = slingshot.1to18$slingPseudotime_1,
y = forcats::fct_relevel(slingshot.1to18$SCT_snn_res.0.8,
"1", "5", "13", "15", "18" , "21"),
colour = SCT_snn_res.0.8)) +
geom_quasirandom(groupOnX = F) +
facet_grid(rows = vars(genotype)) +
scale_color_tableau() + theme_classic() +
xlab("Slingshot pseudotime") + ylab("Clusters") +
ggtitle("") + NoLegend()## Warning: Removed 30 rows containing missing values (position_quasirandom).## Warning: Removed 81 rows containing missing values (position_quasirandom).grid.text("(Fig22)", 0.15, 0.85)
\(~\)
\(~\)
Statistical test between distribution of pseudotime1 and 2 in wild-type and wnt11
# On pseudotime one
ks.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])## Warning in ks.test(slingPseudotime(slingshot.wild.1to18)[, 1],
## slingPseudotime(slingshot.wnt11.1to18)[, : p-value will be approximate in the
## presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## D = 0.10589, p-value = 0.0001303
## alternative hypothesis: two-sided# On pseudotime two
ks.test(slingPseudotime(slingshot.wild.1to18)[, 2], slingPseudotime(slingshot.wnt11.1to18)[, 2])## Warning in ks.test(slingPseudotime(slingshot.wild.1to18)[, 2],
## slingPseudotime(slingshot.wnt11.1to18)[, : p-value will be approximate in the
## presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: slingPseudotime(slingshot.wild.1to18)[, 2] and slingPseudotime(slingshot.wnt11.1to18)[, 2]
## D = 0.099284, p-value = 0.000381
## alternative hypothesis: two-sided\(~\)
\(~\)
Boxplots
summary(slingPseudotime(slingshot.wild.1to18))[,1]##
## "Min. : 0.000 " "1st Qu.: 1.019 " "Median : 2.841 " "Mean : 3.275 "
##
## "3rd Qu.: 5.228 " "Max. :10.365 " "NA's :35 "summary(slingPseudotime(slingshot.wnt11.1to18))[,1]##
## "Min. : 0.000 " "1st Qu.: 1.525 " "Median : 3.376 " "Mean : 3.949 "
##
## "3rd Qu.: 6.012 " "Max. :10.199 " "NA's :80 "boxplot(list(wild = slingPseudotime(slingshot.wild.1to18)[,1] , wnt11 = slingPseudotime(slingshot.wnt11.1to18)[,1]) ,
horizontal = T,
notch = T, # Confidence interval for median
main = "Pseudotime One",
sub = "",
xlab = "",
col = "gray" )
grid.text("(Fig23)", 0.15, 0.85)
boxplot(list(wild = slingPseudotime(slingshot.wild.1to18)[,2] , wnt11 = slingPseudotime(slingshot.wnt11.1to18)[,2]) ,
horizontal = T,
notch = T, # Confidence interval for median
main = "Pseudotime two",
sub = "",
xlab = "",
col = "gray50" )
grid.text("(Fig24)", 0.15, 0.85)
\(~\)
\(~\)
QQplots
QQplot shows that there is a inharmony between the wild-type and wnt11 cells in the later pseudotimes.
par(mfrow = c(1,2))
qqplot(slingPseudotime(slingshot.wild.1to18)[,1] , slingPseudotime(slingshot.wnt11.1to18)[,1] ,
main = "Quantile-Quantile Plot Psudotime One",
xlab = "wild-type Percentiles",
ylab = "Wnt11 Percentiles",
col.axis = 'blue', col.lab = 'red', cex.axis = 1, cex.lab = 1 , cex.main = 1)
grid.text("(Fig25)", 0.15, 0.85)
qqplot(slingPseudotime(slingshot.wild.1to18)[,2] , slingPseudotime(slingshot.wnt11.1to18)[,2] ,
main = "Quantile-Quantile Plot Psudotime two",
xlab = "wild-type Percentiles",
ylab = "Wnt11 Percentiles",
col.axis = 'blue', col.lab = 'red', cex.axis = 1, cex.lab = 1 , cex.main = 1)
grid.text("(Fig26)", 0.25, 0.85)
dev.off()## null device
## 1\(~\)
\(~\)
Density Plots
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,1])
plot(dens, frame = F , col = "steelblue", main = "Wild-type Psudotime One")
polygon(dens, col = "steelblue")
grid.text("(Fig27)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,1])
plot(dens, frame = FALSE, col = "indianred1", main = "Wnt11 Psudotime One")
polygon(dens, col = "indianred1")
grid.text("(Fig28)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,1])
plot(dens, frame = F , col = "steelblue", main = "Psudotime One" , xlab = "")
grid.text("(Fig29)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,1])
lines(dens, col = "indianred1", main = "")
\(~\)
\(~\)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,2])
plot(dens, frame = F , col = "steelblue", main = "Wild-type Psudotime Two")
polygon(dens, col = "steelblue")
grid.text("(Fig30)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,2])
plot(dens, frame = FALSE, col = "indianred1", main = "Wnt11 Psudotime Two")
polygon(dens, col = "indianred1")
grid.text("(Fig31)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wild.1to18)))[,2])
plot(dens, frame = F , col = "steelblue", main = "Psudotime Two" , xlab = "")
grid.text("(Fig32)", 0.15, 0.75)
dens = density(as.data.frame(na.omit(slingPseudotime(slingshot.wnt11.1to18)))[,2])
lines(dens, col = "indianred1", main = "")
\(~\)
\(~\)
Differential Expression Analysis
For each condition (wild and wnt11), a smooth average expression profile along pseudotime will be estimated for each gene, using a negative binomial generalized additive model (NB-GAM). But first, we need to Evaluate the optimal number of knots required for fitGAM.
plan("multiprocess", workers = 18)## Warning in supportsMulticoreAndRStudio(...): [ONE-TIME WARNING] Forked
## processing ('multicore') is not supported when running R from RStudio
## because it is considered unstable. For more details, how to control forked
## processing or not, and how to silence this warning in future R sessions, see ?
## parallelly::supportsMulticoreBPPARAMS <- BiocParallel::bpparam()
BPPARAMS$workers <- 12icMat = evaluateK(slingshot.1to18, conditions = factor(colData(slingshot.1to18)$genotype,
ordered = FALSE), nGenes = 500, k = 3:10, parallel=TRUE, BPPARAM = BPPARAMS)
Figure33
\(~\)
In the left panel, boxplots, the skewness in the distribution of gene-level AIC values changes as the number of knots increases from eight knots. In addition, the Average AIC levels off as it reaches the eighth knot.
\(~\)
\(~\)
GAM.1to18 <- fitGAM(counts = slingshot.1to18@assays@data$counts %>% as.matrix(), sds = sds.1to18, conditions = factor(colData(slingshot.1to18)$genotype, ordered = FALSE), nknots = 8)\(~\)
\(~\)
Perform statistical test to check whether gene expression is constant across pseudotime within a lineage
To assess significant changes in gene expression as a function of pseudotime within each lineage, we use the “associationTest” (Wald test), which tests the null hypothesis that gene expression is not a function of pseudotime, i.e., whether the estimated smoothers are significantly varying as a function of pseudotime within each lineage. The lineages=TRUE argument specifies that we would like the results for each lineage separately, asides from the default global test, which tests for significant associations across all lineages in the trajectory simultaneously. Further, we specify a log2 fold change cut-off to test against using the l2fc argument.
rowData(GAM.1to18)$assocRes <- associationTest(GAM.1to18, lineages = TRUE, l2fc = log2(2))assocRes <- rowData(GAM.1to18)$assocRes
all(assocRes$meanLogFC > 0)## [1] TRUE\(~\)
\(~\)
Selecting genes with fdr-adjusted p-value less than five percent on lineage1 for both conditions.
wildGenes <- rownames(assocRes)[
which(p.adjust(assocRes$`pvalue_lineage1_conditionwild-type`, "fdr") <= 0.05)
]
wnt11Genes <- rownames(assocRes)[
which(p.adjust(assocRes$`pvalue_lineage1_conditionwnt11-hom`, "fdr") <= 0.05)
]
length(wildGenes) ## [1] 391length(wnt11Genes) ## [1] 358\(~\)
\(~\)
Genes that are common between wild and wnt11
intersect(wildGenes , wnt11Genes)## [1] "Gsta3" "Col3a1" "Col5a2" "Fn1"
## [5] "Ugt1a9" "Arl4c" "Tnfrsf11a" "Tmem37"
## [9] "Plekha6" "Cfh" "4930500M09Rik" "Cenpf"
## [13] "Irf6" "Vim" "Pax8" "Gsn"
## [17] "Zeb2" "Dpp4" "Ttn" "Serping1"
## [21] "Rapsn" "Gm10804" "Rcn1" "Prnp"
## [25] "Myl9" "Mafb" "Ube2c" "Bgn"
## [29] "Cited1" "Maged2" "Tmsb4x" "Ccna2"
## [33] "Gucy1a1" "Fga" "S100a6" "S100a10"
## [37] "Ecm1" "Hist2h2ac" "Pdzk1" "Olfml3"
## [41] "Arhgap29" "Abca4" "Npnt" "Penk"
## [45] "Tpm2" "Gm12678" "Slc5a9" "Cdc20"
## [49] "Cdca8" "Eva1b" "Gm12999" "Sfn"
## [53] "Aunip" "Mfap2" "Chd5" "Mmp23"
## [57] "Ttll10" "Cenpa" "Emilin1" "Ldb2"
## [61] "Cckar" "Pdgfra" "Igfbp7" "Gm33050"
## [65] "Bmp3" "Plac8" "Arhgap24" "Pcolce"
## [69] "Rasl11a" "Col1a2" "Cald1" "Tmem140"
## [73] "Rarres2" "Npy" "Gprin3" "Ndnf"
## [77] "Asprv1" "Bcam" "Cblc" "Prodh2"
## [81] "Fxyd7" "Klk10" "Emp3" "Fah"
## [85] "AI314278" "Prss23" "Stard10" "Rassf10"
## [89] "Acsm1" "Crym" "Tgfb1i1" "H19"
## [93] "Cdkn1c" "Ccnd1" "Akap12" "Gm47889"
## [97] "Upb1" "Cnn2" "Timp3" "Srgap1"
## [101] "Eef1akmt3" "Hmgb2" "BC030500" "Sh2d4a"
## [105] "Pde4c" "Elmo3" "Clec18a" "Snai3"
## [109] "D830030K20Rik" "Bmp4" "Clu" "Pdlim2"
## [113] "Izumo1r" "Ntm" "Tagln" "Hspb2"
## [117] "Cryab" "Slc35f2" "Islr2" "Pclaf"
## [121] "Tpm1" "Anxa2" "Dag1" "Selenom"
## [125] "4921536K21Rik" "Emid1" "Col23a1" "Pmp22"
## [129] "Tmem102" "Pimreg" "Hic1" "Lhx1"
## [133] "Lhx1os" "Car4" "Col1a1" "Top2a"
## [137] "Kcnj16" "Birc5" "Hist1h1b" "Hist1h2ap"
## [141] "Hist1h2ae" "Hist1h4d" "Hist1h1e" "Hist1h2ab"
## [145] "Hist1h1a" "Gmnn" "Gm17750" "Bhmt2"
## [149] "Ccnb1" "Rrm2" "Scin" "Mis18bp1"
## [153] "Meg3" "Agxt2" "Npr3" "Osr2"
## [157] "Cthrc1" "Lgals1" "Igfbp6" "Cldn5"
## [161] "Stfa3" "Fstl1" "Pcp4" "Ripk4"
## [165] "Clic5" "Uhrf1" "Lbh" "Cdc42ep3"
## [169] "Six2" "Dsg2" "Mep1b" "Ccdc178"
## [173] "Cd248" "Slc22a19" "Cpn1" "Emx2"\(~\)
\(~\)
Visualization of number of DE genes that present in either both or one cell type.
UpSetR::upset(fromList(list(wild = wildGenes, wnt11 = wnt11Genes)))
grid.text("(Fig34)", 0.15, 0.85)
\(~\)
\(~\)
DEGs that present only in wild-type
wildGenes[!(wildGenes %in% wnt11Genes)]## [1] "Gm16894" "8430432A02Rik" "Plcd4" "Inpp5d"
## [5] "Gm26720" "Cdh20" "Rassf5" "Gm29629"
## [9] "Fmod" "Mybph" "Myog" "Ppfia4"
## [13] "Rgs16" "Pbx1" "Slamf1" "Tlr5"
## [17] "Prox1os" "Cd46" "Acbd7" "Cacna1b"
## [21] "Tor4a" "Tubb4b" "Gm996" "Lcn2"
## [25] "Ifih1" "Spc25" "Rapgef4" "Fibin"
## [29] "Rhov" "Nusap1" "Gm14133" "Tpx2"
## [33] "Hnf4aos" "Aurka" "Gpc3" "Xlr4b"
## [37] "Trpc5" "Lhfpl1" "Frmpd4" "Fabp4"
## [41] "Gm5103" "Mme" "Golim4" "Gm15998"
## [45] "Hfe2" "5330417C22Rik" "Fam166b" "Stra6l"
## [49] "Smc2" "Adamtsl1" "Kif2c" "Ccdc24"
## [53] "Zc3h12a" "Matn1" "Ncmap" "Sh2d5"
## [57] "2700016F22Rik" "Slc25a34" "Gpr157" "Gm26840"
## [61] "Tmem88b" "Lrrd1" "Gm8773" "Sema3a"
## [65] "Fgfbp1" "Nipal1" "Ptpn13" "Abcg3"
## [69] "Idua" "Plcxd1" "Selplg" "Mlxipl"
## [73] "Fzd9" "Gm15498" "Gper1" "Mmd2"
## [77] "2900089D17Rik" "Irf5" "Rncr4" "Trh"
## [81] "Frmd4b" "Syn2" "Cxcl12" "D330020A13Rik"
## [85] "Gm38910" "Apold1" "Bhlhe41" "Pnmal1"
## [89] "Cd177" "Rinl" "Gm26810" "Nphs1"
## [93] "Cd33" "Gm45669" "Kcnj14" "1700011D18Rik"
## [97] "Prc1" "Cemip" "Me3" "Hbb-bt"
## [101] "Hbb-bs" "Gm26690" "D7Ertd443e" "Gm44732"
## [105] "Rgs17" "Adgrg6" "Cenpw" "Mfsd4b3"
## [109] "Tspan15" "Cdk1" "Plk5" "Gna15"
## [113] "Gm32687" "Eid3" "Ntn4" "Myl6"
## [117] "Efnb2" "Proscos" "Mfap3l" "Gm15991"
## [121] "Lzts1" "Iqcm" "Gm10645" "Gm45767"
## [125] "Gm3636" "Gm3739" "Duxbl3" "Stab1"
## [129] "Eddm3b" "Rnase6" "Trac" "Cma1"
## [133] "Loxl2" "Itm2b" "Igsf9b" "Gm15520"
## [137] "Adamts15" "Slc37a2" "Gm10658" "9230112J17Rik"
## [141] "Unc13c" "Gcm1" "Fam46a" "Fbxw15"
## [145] "Dlec1" "Ccdc13" "Cdcp1" "Plek"
## [149] "Hba-a1" "Hba-a2" "Hmmr" "Gabrb2"
## [153] "Sox30" "Flt4" "Hs3st3b1" "Aurkb"
## [157] "Tm4sf5" "Rilp" "Ccl12" "Wfdc17"
## [161] "Sept4" "Dgkeos" "Sp6" "Csf3"
## [165] "Gast" "Icam2" "Pecam1" "Abca6"
## [169] "Fn3krp" "Akr1c21" "Hist1h3h" "Hist1h2bc"
## [173] "Agtr1a" "Irf4" "Gm26514" "A330048O09Rik"
## [177] "Pfn3" "Gm48488" "Ctsl" "Gm47123"
## [181] "BC048507" "Dhfr" "Tmem171" "Gm10734"
## [185] "Fam228a" "Srp54c" "Tdrd9" "A530016L24Rik"
## [189] "Rgs22" "Pick1.1" "Gm16059" "1810041L15Rik"
## [193] "Gm26798" "Shank3" "Pde1b" "Ciita"
## [197] "Aifm3" "Rtl10" "Il1rap" "Epha3"
## [201] "Cldn6" "Prss27" "Slc9a3r2" "Itpr3"
## [205] "Scube3" "Platr17" "Trem2" "Apobec2"
## [209] "Arhgap28" "Mkx" "Rnf125" "Cd14"
## [213] "Gnal" "Rin1" "Trpm3"\(~\)
DEGs that present only in wnt11
wnt11Genes[!(wnt11Genes %in% wildGenes)]## [1] "Eya1" "D630023F18Rik" "Pid1" "Col6a3"
## [5] "Cdk18" "Tmem81" "Angptl1" "Atp1b1"
## [9] "Gm37422" "Rd3" "Nrarp" "Phyhd1"
## [13] "Pip5kl1" "Lzts3" "Plcb1" "Jag1"
## [17] "Hnf4a" "5031425F14Rik" "Cybb" "Rtl8a"
## [21] "Rtl8b" "Gm14827" "Gpm6b" "Cldn11"
## [25] "Anxa5" "Tm4sf1" "Enpep" "Lef1"
## [29] "Arid3c" "Col15a1" "Acnat1" "Aldob"
## [33] "Kif12" "Nfib" "Ccdc17" "Kcnq4"
## [37] "Tmem54" "Gm12976" "Fam229a" "Wnt4"
## [41] "Hes5" "Tmem52" "Gm26894" "Gm43660"
## [45] "Tacc3" "Gm15477" "Gabra2" "Afm"
## [49] "Cdkl2" "Tmem150c" "Gm26711" "Oasl2"
## [53] "Wnt16" "Lrrc4" "Podxl" "Svopl"
## [57] "Epha1" "Aoc1" "Tacstd2" "Eva1a"
## [61] "Mitf" "Pparg" "Clec4a2" "Slco1a6"
## [65] "Cdc42ep5" "Phldb3" "Aspdh" "Slc17a7"
## [69] "Chrna7" "Agbl1" "Wdr93" "Gm44949"
## [73] "Tm6sf1" "Slco2b1" "Hbb-y" "Spon1"
## [77] "Atp2a1" "Aldoa.1" "Gm44759" "Ifitm3"
## [81] "Gm16982" "Cd81" "Lama2" "Lama4"
## [85] "Lilrb4a" "Srgn" "Col6a2" "Dcn"
## [89] "Atp2b1" "Alx1" "Lrriq1" "Nav3"
## [93] "Ptprr" "Gm36041" "Mcemp1" "Tgfbr3l"
## [97] "Col4a1" "Proz" "Gm32050" "Enpp6"
## [101] "Lrrc25" "Il15" "Nfix" "Cdh16"
## [105] "Gm33023" "Slc35f3" "Cfap70" "Gm30108"
## [109] "Wnt5a" "Vstm4" "Arhgap22" "Nrg3"
## [113] "Fermt2" "Gm34934" "1700011H14Rik" "Cebpe"
## [117] "Myh7" "Dct" "Gria4" "Nxpe2"
## [121] "Arhgap20os" "Crabp1" "Islr" "Trpc1"
## [125] "Als2cl" "Hba-x" "Ebf1" "Acsl6"
## [129] "Anxa6" "Slc5a10" "Tnk1" "Slc6a4"
## [133] "Hnf1b" "Hoxb1" "Nags" "Cygb"
## [137] "Gm11728" "CT010465.1" "Nid1" "Hist1h1d"
## [141] "Hist1h3e" "Hist1h3c" "Dcdc2a" "Dsp"
## [145] "Rnf144b" "Id4" "Cks2" "Cdhr2"
## [149] "Gas1" "Zfp455" "Foxd1" "Gpx8"
## [153] "Arl4a" "Rhoj" "Gm26571" "Fbln5"
## [157] "AC163040.1" "Dlk1" "Amn" "Rai14"
## [161] "Pdzd2" "Snhg18" "Ncald" "Galr3"
## [165] "B130046B21Rik" "Smagp" "Cdc45" "Nrros"
## [169] "Hcls1" "Sytl3" "Dll1" "Tmem8"
## [173] "Grm4" "Gm16279" "Ttbk1" "1700061G19Rik"
## [177] "Pdgfrb" "Ablim3" "Myo5b" "Acta2"
## [181] "Ifit3" "Scd3"\(~\)
\(~\)
Heatmap of wild DE genes based on mean smoother. It Gets smoothers estimated by tradeSeq along a grid. This function does not return fitted values but rather the predicted mean smoother, for a user-defined grid of points. If tidy = FALSE and the trajectory consists of 2 lineages and nPoints=100, then the returned matrix will have 2*100 columns, where the first 100 correspond to the first lineage and columns 101-200 to the second lineage.
yhatSmooth <- predictSmooth(GAM.1to18, gene = wildGenes, nPoints = 100, tidy = FALSE)
heatSmooth <- pheatmap(t(scale(t(yhatSmooth[, 1:100]))),
cluster_cols = F, show_rownames = T,
show_colnames = F, fontsize_row = 2 , main = "Wild-type DEGs on lineage1")
grid.text("(Fig35)", 0.15, 0.85)
\(~\)
\(~\)
Dendrogram of predicted smooth average expression
c1 <- sort(cutree(heatSmooth$tree_row, k = 6))
c1## Gsta3 Inpp5d Arl4c Tnfrsf11a Tmem37
## 1 1 1 1 1
## Rassf5 Plekha6 Fmod Mybph Myog
## 1 1 1 1 1
## Rgs16 Slamf1 Tlr5 Irf6 Acbd7
## 1 1 1 1 1
## Gsn Dpp4 Ifih1 Ttn Rapsn
## 1 1 1 1 1
## Rcn1 Rhov Prnp Gm14133 Xlr4b
## 1 1 1 1 1
## Fabp4 Mme 5330417C22Rik Abca4 Npnt
## 1 1 1 1 1
## Fam166b Stra6l Slc5a9 Ccdc24 Gm12999
## 1 1 1 1 1
## Sfn Ncmap 2700016F22Rik Ttll10 Gm8773
## 1 1 1 1 1
## Fgfbp1 Cckar Nipal1 Igfbp7 Arhgap24
## 1 1 1 1 1
## Ptpn13 Abcg3 Idua Plcxd1 Mlxipl
## 1 1 1 1 1
## Gm15498 Gper1 Irf5 Rncr4 Tmem140
## 1 1 1 1 1
## Frmd4b Syn2 D330020A13Rik Apold1 Cblc
## 1 1 1 1 1
## Gm26810 Klk10 Fah Prss23 Me3
## 1 1 1 1 1
## Stard10 Hbb-bt Hbb-bs Gm26690 D7Ertd443e
## 1 1 1 1 1
## Gm44732 Rgs17 Adgrg6 Gm47889 Upb1
## 1 1 1 1 1
## Plk5 Eef1akmt3 Proscos BC030500 Mfap3l
## 1 1 1 1 1
## Gm15991 Pde4c Iqcm Elmo3 Clec18a
## 1 1 1 1 1
## Snai3 D830030K20Rik Gm3636 Gm3739 Stab1
## 1 1 1 1 1
## Bmp4 Eddm3b Rnase6 Trac Cma1
## 1 1 1 1 1
## Clu Itm2b Izumo1r Adamts15 Slc37a2
## 1 1 1 1 1
## Hspb2 Slc35f2 Gm10658 Unc13c Dag1
## 1 1 1 1 1
## Fbxw15 Dlec1 Cdcp1 Emid1 Hba-a1
## 1 1 1 1 1
## Hba-a2 Gabrb2 Hs3st3b1 Tm4sf5 Car4
## 1 1 1 1 1
## Sept4 Sp6 Abca6 Kcnj16 Fn3krp
## 1 1 1 1 1
## Hist1h3h Agtr1a A330048O09Rik Ctsl BC048507
## 1 1 1 1 1
## Gm17750 Tmem171 Gm10734 Scin Osr2
## 1 1 1 1 1
## Rgs22 Cthrc1 1810041L15Rik Pde1b Aifm3
## 1 1 1 1 1
## Il1rap Ripk4 Cldn6 Prss27 Itpr3
## 1 1 1 1 1
## Platr17 Trem2 Apobec2 Mkx Dsg2
## 1 1 1 1 1
## Ccdc178 Cd14 Trpm3 Cpn1 Emx2
## 1 1 1 1 1
## Gm16894 8430432A02Rik Plcd4 Gm26720 Cdh20
## 2 2 2 2 2
## Ppfia4 4930500M09Rik Prox1os Cd46 Cacna1b
## 2 2 2 2 2
## Tor4a Gm996 Lcn2 Rapgef4 Cited1
## 2 2 2 2 2
## Trpc5 Lhfpl1 Frmpd4 Gm5103 Gm15998
## 2 2 2 2 2
## Hfe2 Tpm2 Zc3h12a Matn1 Sh2d5
## 2 2 2 2 2
## Mfap2 Slc25a34 Gpr157 Chd5 Gm26840
## 2 2 2 2 2
## Tmem88b Lrrd1 Selplg Fzd9 Mmd2
## 2 2 2 2 2
## 2900089D17Rik Gprin3 Asprv1 Bhlhe41 Pnmal1
## 2 2 2 2 2
## Cd177 Rinl Gm45669 Kcnj14 Emp3
## 2 2 2 2 2
## 1700011D18Rik Cemip Crym Gna15 Gm32687
## 2 2 2 2 2
## Eid3 Lzts1 Gm10645 Duxbl3 Igsf9b
## 2 2 2 2 2
## Islr2 9230112J17Rik Fam46a Ccdc13 4921536K21Rik
## 2 2 2 2 2
## Plek Sox30 Flt4 Rilp Ccl12
## 2 2 2 2 2
## Wfdc17 Dgkeos Csf3 Gast Icam2
## 2 2 2 2 2
## Irf4 Gm26514 Pfn3 Gm48488 Gm47123
## 2 2 2 2 2
## Fam228a Srp54c Tdrd9 A530016L24Rik Gm16059
## 2 2 2 2 2
## Gm26798 Shank3 Ciita Stfa3 Six2
## 2 2 2 2 2
## Rnf125 Gnal Rin1 Col3a1 Col5a2
## 2 2 2 3 3
## Fn1 Cfh Pbx1 Zeb2 Serping1
## 3 3 3 3 3
## Fibin Myl9 Gpc3 Bgn Maged2
## 3 3 3 3 3
## Tmsb4x Gucy1a1 S100a6 S100a10 Ecm1
## 3 3 3 3 3
## Olfml3 Penk Adamtsl1 Eva1b Sema3a
## 3 3 3 3 3
## Emilin1 Ldb2 Pdgfra Plac8 Pcolce
## 3 3 3 3 3
## Col1a2 Cald1 Rarres2 Fxyd7 Tgfb1i1
## 3 3 3 3 3
## Cdkn1c Cnn2 Pdlim2 Ntm Anxa2
## 3 3 3 3 3
## Selenom Col23a1 Hic1 Col1a1 Meg3
## 3 3 3 3 3
## Lgals1 Epha3 Lbh Cdc42ep3 Cd248
## 3 3 3 3 3
## Ugt1a9 Gm29629 Gm10804 Hnf4aos Fga
## 4 4 4 4 4
## Pdzk1 Gm12678 Gm33050 Trh Gm38910
## 4 4 4 4 4
## Prodh2 Cd33 AI314278 Acsm1 Mfsd4b3
## 4 4 4 4 4
## Gm45767 Gm15520 Gcm1 Tmem102 Pecam1
## 4 4 4 4 4
## Akr1c21 Hist1h2bc Bhmt2 Agxt2 Pick1.1
## 4 4 4 4 4
## Rtl10 Mep1b Slc22a19 Cenpf Tubb4b
## 4 4 4 5 5
## Spc25 Nusap1 Tpx2 Ube2c Aurka
## 5 5 5 5 5
## Ccna2 Hist2h2ac Smc2 Kif2c Cdc20
## 5 5 5 5 5
## Cdca8 Aunip Cenpa Prc1 Cenpw
## 5 5 5 5 5
## Cdk1 Hmgb2 Pclaf Hmmr Aurkb
## 5 5 5 5 5
## Pimreg Top2a Birc5 Hist1h1b Hist1h2ap
## 5 5 5 5 5
## Hist1h2ae Hist1h4d Hist1h1e Hist1h2ab Hist1h1a
## 5 5 5 5 5
## Gmnn Dhfr Ccnb1 Rrm2 Mis18bp1
## 5 5 5 5 5
## Uhrf1 Vim Pax8 Mafb Golim4
## 5 6 6 6 6
## Arhgap29 Mmp23 Bmp3 Rasl11a Npy
## 6 6 6 6 6
## Ndnf Cxcl12 Bcam Nphs1 Rassf10
## 6 6 6 6 6
## H19 Ccnd1 Akap12 Tspan15 Timp3
## 6 6 6 6 6
## Ntn4 Srgap1 Myl6 Efnb2 Sh2d4a
## 6 6 6 6 6
## Loxl2 Tagln Cryab Tpm1 Pmp22
## 6 6 6 6 6
## Lhx1 Lhx1os Npr3 Igfbp6 Cldn5
## 6 6 6 6 6
## Fstl1 Pcp4 Slc9a3r2 Scube3 Clic5
## 6 6 6 6 6
## Arhgap28
## 6table(c1)## c1
## 1 2 3 4 5 6
## 150 88 47 28 38 40plot(heatSmooth$tree_row , sub = "" , main = "" , xlab = "" , cex = 0.3)
grid.text("(Fig36)", 0.15, 0.85)
\(~\)
\(~\)
Visualizing the heatmap for the actual fitted values in sorted cells based on pseudotimes. In other words, cells with larger pseudotime are toward the right-side of the heatmap. In addition, Clustering were applied only on rows (genes).
conditions <- GAM.1to18$tradeSeq$conditions
pt1 <- slingshot.1to18$slingPseudotime_1
## Based on fitted values
yhatCell <- predictCells(GAM.1to18, gene=wildGenes)
dim(yhatCell)## [1] 391 1847yhatCellwild <- yhatCell[ , conditions == "wild-type"]
yhatCellwild = yhatCellwild[ , !is.na(pt1[conditions == "wild-type"])]
dim(yhatCellwild)## [1] 391 797pt1.wild = pt1[conditions == "wild-type"][!is.na(pt1[conditions == "wild-type"])]
# Order according to pseudotime
# Remove NA pseudotimes
# Separate cells that are in pseudotime zero from the other cells. Then, combines them together to have zero-pseudotime cells at the begining of the matrix, then sorted-non-zero-psuedotime cells.
sum(pt1.wild == 0)## [1] 72yhatCellwildzero = yhatCellwild[ , pt1.wild == 0]
yhatCellwild.nonzero = yhatCellwild[ , pt1.wild != 0]
dim(yhatCellwild.nonzero)## [1] 391 725dim(yhatCellwildzero)## [1] 391 72dim(yhatCellwild)## [1] 391 797pt1.wild.zero = pt1.wild[pt1.wild == 0]
pt1.wild.nonzero = pt1.wild[pt1.wild != 0]
pt1.wild.nonzero.sorted = sort(pt1.wild.nonzero, decreasing=FALSE)
yhatCellwild.nonzero.sorted = yhatCellwild.nonzero[ , match(pt1.wild.nonzero.sorted, pt1.wild.nonzero)]
mat = cbind(yhatCellwildzero , yhatCellwild.nonzero.sorted)
dim(mat)## [1] 391 797pheatmap(t(scale(t(mat))), cluster_cols = FALSE,
show_rownames = T, show_colnames=FALSE , fontsize_row = 2 , main = "Wild-type DEGs across lineage1, wild-type sorted cells based on lineage 1")
grid.text("(Fig37)", 0.15, 0.85)
\(~\)
\(~\)
GSEA based on wild-type pseudotimes on lineage1
## C5 category is according to gene ontology grouping: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4707969/pdf/nihms-743907.pdf
geneSets <- msigdbr(species = "Mus musculus", category = "C5", subcategory = "BP")
### filter background to only include genes that we assessed.
geneSets <- geneSets[geneSets$gene_symbol %in% names(GAM.1to18),]
m_list <- geneSets %>% split(x = .$gene_symbol, f = .$gs_name)
stats <- assocRes$`waldStat_lineage1_conditionwild-type`
names(stats) <- rownames(assocRes)
eaRes <- fgsea(pathways = m_list, stats = stats, nperm = 5e4, minSize = 10)## Warning in fgsea(pathways = m_list, stats = stats, nperm = 50000, minSize = 10):
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (8.43% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.## Warning in fgseaSimple(...): There were 26 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic valueseaRes = eaRes %>% arrange(pval)
head(eaRes[,c(1:3,8)] , 40)## pathway
## 1: GOBP_NEGATIVE_REGULATION_OF_VIRAL_INDUCED_CYTOPLASMIC_PATTERN_RECOGNITION_RECEPTOR_SIGNALING_PATHWAY
## 2: GOBP_RIG_I_SIGNALING_PATHWAY
## 3: GOBP_REGULATION_OF_RIG_I_SIGNALING_PATHWAY
## 4: GOBP_REGULATION_OF_VIRAL_INDUCED_CYTOPLASMIC_PATTERN_RECOGNITION_RECEPTOR_SIGNALING_PATHWAY
## 5: GOBP_NEGATIVE_REGULATION_OF_DEFENSE_RESPONSE_TO_VIRUS
## 6: GOBP_NEGATIVE_REGULATION_OF_TRANSCRIPTION_BY_COMPETITIVE_PROMOTER_BINDING
## 7: GOBP_HISTONE_METHYLATION
## 8: GOBP_MYD88_INDEPENDENT_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY
## 9: GOBP_RESPONSE_TO_LEUKEMIA_INHIBITORY_FACTOR
## 10: GOBP_POSITIVE_REGULATION_OF_RRNA_PROCESSING
## 11: GOBP_NON_MOTILE_CILIUM_ASSEMBLY
## 12: GOBP_NEGATIVE_REGULATION_OF_MYOTUBE_DIFFERENTIATION
## 13: GOBP_POSITIVE_REGULATION_OF_TRANSCRIPTION_ELONGATION_FROM_RNA_POLYMERASE_II_PROMOTER
## 14: GOBP_TOLL_LIKE_RECEPTOR_4_SIGNALING_PATHWAY
## 15: GOBP_NEGATIVE_REGULATION_OF_DEFENSE_RESPONSE
## 16: GOBP_CYTOPLASMIC_MICROTUBULE_ORGANIZATION
## 17: GOBP_CEREBRAL_CORTEX_NEURON_DIFFERENTIATION
## 18: GOBP_HISTONE_H3_K27_METHYLATION
## 19: GOBP_NEGATIVE_REGULATION_OF_RESPONSE_TO_BIOTIC_STIMULUS
## 20: GOBP_INOSITOL_TRISPHOSPHATE_METABOLIC_PROCESS
## 21: GOBP_LIPOPOLYSACCHARIDE_MEDIATED_SIGNALING_PATHWAY
## 22: GOBP_REGULATION_OF_RRNA_PROCESSING
## 23: GOBP_POSITIVE_REGULATION_OF_SMOOTH_MUSCLE_CELL_DIFFERENTIATION
## 24: GOBP_POSITIVE_REGULATION_OF_NEUROTRANSMITTER_SECRETION
## 25: GOBP_POSITIVE_REGULATION_OF_INTERFERON_GAMMA_PRODUCTION
## 26: GOBP_NUCLEAR_MEMBRANE_ORGANIZATION
## 27: GOBP_CIRCADIAN_REGULATION_OF_GENE_EXPRESSION
## 28: GOBP_RESPONSE_TO_MAGNESIUM_ION
## 29: GOBP_REGULATION_OF_NUCLEAR_TRANSCRIBED_MRNA_CATABOLIC_PROCESS_DEADENYLATION_DEPENDENT_DECAY
## 30: GOBP_REGULATION_OF_VASCULAR_ASSOCIATED_SMOOTH_MUSCLE_CELL_DIFFERENTIATION
## 31: GOBP_POSITIVE_REGULATION_OF_G_PROTEIN_COUPLED_RECEPTOR_SIGNALING_PATHWAY
## 32: GOBP_POSITIVE_REGULATION_OF_SMOOTH_MUSCLE_CONTRACTION
## 33: GOBP_REGULATION_OF_LONG_TERM_SYNAPTIC_DEPRESSION
## 34: GOBP_MYD88_DEPENDENT_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY
## 35: GOBP_POSITIVE_REGULATION_OF_SIGNAL_TRANSDUCTION_BY_P53_CLASS_MEDIATOR
## 36: GOBP_POSITIVE_REGULATION_OF_NEUROTRANSMITTER_TRANSPORT
## 37: GOBP_POSITIVE_REGULATION_OF_RELEASE_OF_CYTOCHROME_C_FROM_MITOCHONDRIA
## 38: GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_APOPTOTIC_PROCESS
## 39: GOBP_POSITIVE_REGULATION_OF_ALCOHOL_BIOSYNTHETIC_PROCESS
## 40: GOBP_POSITIVE_REGULATION_OF_CELLULAR_CARBOHYDRATE_METABOLIC_PROCESS
## pathway
## pval padj leadingEdge
## 1: 0.001138057 0.8789942 Rnf125,Tkfc
## 2: 0.001641642 0.8789942 Rnf125,Ddx60
## 3: 0.001749201 0.8789942 Rnf125,Ddx60
## 4: 0.001804258 0.8789942 Rnf125,Ddx60
## 5: 0.002103703 0.8789942 Rnf125,Mill2
## 6: 0.003080147 0.8789942 Bhlhe41,Muc1
## 7: 0.003379932 0.8789942 Chd5,Smyd1
## 8: 0.003560071 0.8789942 Cd14,Cd40
## 9: 0.003639927 0.8789942 Rnf125,Tex14
## 10: 0.003952569 0.8789942 Trmt112,Wdr43,Riok2,Dimt1,Utp15,Riok1,...
## 11: 0.004399912 0.8789942 Ccdc13,Cep350
## 12: 0.004756417 0.8789942 Bhlhe41,Notch1
## 13: 0.005376344 0.8789942 Tcerg1,Gtf2f1,Brd4,Eapp,Cdk13,Cdk12,...
## 14: 0.005500110 0.8789942 Cd14,Lbp
## 15: 0.005719886 0.8789942 Rnf125,Gper1,Serping1
## 16: 0.006439871 0.8789942 Ccdc13,Mapt
## 17: 0.006504629 0.8789942 Chd5,Emx1
## 18: 0.006674917 0.8789942 Chd5,Phf19
## 19: 0.006779864 0.8789942 Rnf125,Serping1
## 20: 0.008088004 0.8789942 Gper1,P2ry6
## 21: 0.008119838 0.8789942 Cd14,Sigirr
## 22: 0.008415147 0.8789942 Trmt112,Wdr43,Kat2b,Riok2,Dimt1,Utp15,...
## 23: 0.008423258 0.8789942 Gper1,Eng
## 24: 0.010326848 0.8789942 Gper1,Nlgn1
## 25: 0.010379792 0.8789942 Cd14,Wnt5a
## 26: 0.010495205 0.8789942 Gper1,Reep4
## 27: 0.010919782 0.8789942 Bhlhe41,Id4
## 28: 0.011137163 0.8789942 Cd14,Ccnd1
## 29: 0.011220196 0.8789942 Nanos1,Cpeb3,Zfp36l2,Tnrc6b,Ago2,Pabpc1,...
## 30: 0.011320984 0.8789942 Gper1,Eng
## 31: 0.011399964 0.8789942 Gper1,Tmod2
## 32: 0.011525642 0.8789942 Gper1,Oxtr
## 33: 0.012215820 0.8789942 Shank3,Ager
## 34: 0.012280491 0.8789942 Cd14,Tlr5,Cd36
## 35: 0.012481968 0.8789942 Chd5,Hic1
## 36: 0.013531174 0.8789942 Gper1,Itgb1
## 37: 0.013602448 0.8789942 Gper1,Hrk
## 38: 0.014905591 0.8789942 Gper1,Cd248
## 39: 0.015216866 0.8789942 Gper1,Wnt4
## 40: 0.015219696 0.8789942 Gper1,Actn3
## pval padj leadingEdge\(~\)
\(~\)
GSEA based on wnt11 pseudotimes on lineage1
stats <- assocRes$`waldStat_lineage1_conditionwnt11-hom`
names(stats) <- rownames(assocRes)
eaRes <- fgsea(pathways = m_list, stats = stats, nperm = 5e4, minSize = 10)## Warning in fgsea(pathways = m_list, stats = stats, nperm = 50000, minSize = 10):
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (5.94% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.## Warning in fgseaSimple(...): There were 14 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic valueseaRes = eaRes %>% arrange(pval)
head(eaRes[,c(1:3,8)] , 40)## pathway
## 1: GOBP_VESICLE_MEDIATED_TRANSPORT_IN_SYNAPSE
## 2: GOBP_SYNAPTIC_VESICLE_RECYCLING
## 3: GOBP_PRESYNAPTIC_ENDOCYTOSIS
## 4: GOBP_REGULATION_OF_SYNAPTIC_VESICLE_RECYCLING
## 5: GOBP_REGULATION_OF_SYNAPTIC_VESICLE_ENDOCYTOSIS
## 6: GOBP_PHOSPHATE_ION_TRANSPORT
## 7: GOBP_FATTY_ACID_TRANSMEMBRANE_TRANSPORT
## 8: GOBP_L_GLUTAMATE_TRANSMEMBRANE_TRANSPORT
## 9: GOBP_POSITIVE_REGULATION_OF_RRNA_PROCESSING
## 10: GOBP_ACIDIC_AMINO_ACID_TRANSPORT
## 11: GOBP_GLUTAMATE_SECRETION
## 12: GOBP_L_ALPHA_AMINO_ACID_TRANSMEMBRANE_TRANSPORT
## 13: GOBP_L_AMINO_ACID_TRANSPORT
## 14: GOBP_AMINO_ACID_TRANSMEMBRANE_TRANSPORT
## 15: GOBP_AMINO_ACID_TRANSPORT
## 16: GOBP_ORGANIC_ACID_TRANSMEMBRANE_TRANSPORT
## 17: GOBP_POSITIVE_REGULATION_OF_TRANSCRIPTION_ELONGATION_FROM_RNA_POLYMERASE_II_PROMOTER
## 18: GOBP_REGULATION_OF_RRNA_PROCESSING
## 19: GOBP_CELLULAR_RESPONSE_TO_ETHANOL
## 20: GOBP_DICARBOXYLIC_ACID_TRANSPORT
## 21: GOBP_MAMMARY_GLAND_EPITHELIAL_CELL_DIFFERENTIATION
## 22: GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_PEPTIDE_ANTIGEN_VIA_MHC_CLASS_I
## 23: GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION_OF_EXOGENOUS_PEPTIDE_ANTIGEN_VIA_MHC_CLASS_I
## 24: GOBP_CELL_REDOX_HOMEOSTASIS
## 25: GOBP_POSITIVE_REGULATION_OF_EXCITATORY_POSTSYNAPTIC_POTENTIAL
## 26: GOBP_NEUROTRANSMITTER_TRANSPORT
## 27: GOBP_SYNAPTIC_TRANSMISSION_DOPAMINERGIC
## 28: GOBP_C_TERMINAL_PROTEIN_AMINO_ACID_MODIFICATION
## 29: GOBP_UTP_METABOLIC_PROCESS
## 30: GOBP_G_PROTEIN_COUPLED_GLUTAMATE_RECEPTOR_SIGNALING_PATHWAY
## 31: GOBP_RESPIRATORY_BURST
## 32: GOBP_NEGATIVE_REGULATION_OF_EPIDERMAL_GROWTH_FACTOR_ACTIVATED_RECEPTOR_ACTIVITY
## 33: GOBP_NEUROTRANSMITTER_BIOSYNTHETIC_PROCESS
## 34: GOBP_CALCIUM_ION_REGULATED_EXOCYTOSIS
## 35: GOBP_ECTODERMAL_PLACODE_FORMATION
## 36: GOBP_MODULATION_OF_EXCITATORY_POSTSYNAPTIC_POTENTIAL
## 37: GOBP_INORGANIC_ANION_TRANSPORT
## 38: GOBP_CAMP_METABOLIC_PROCESS
## 39: GOBP_PROTEIN_AUTOPHOSPHORYLATION
## 40: GOBP_SYNAPTIC_TRANSMISSION_GLUTAMATERGIC
## pathway
## pval padj leadingEdge
## 1: 0.0004999900 0.8412843 Slc17a7,Rab11a
## 2: 0.0005399892 0.8412843 Slc17a7,Sncg
## 3: 0.0008199836 0.8412843 Slc17a7,Sncg
## 4: 0.0008608781 0.8412843 Slc17a7,Lrrk2
## 5: 0.0008681257 0.8412843 Slc17a7,Lrrk2
## 6: 0.0020611955 0.8412843 Slc17a7,Slc17a1
## 7: 0.0023400000 0.8412843 Slc17a7,Slc17a8
## 8: 0.0029441807 0.8412843 Slc17a7,Slc17a8
## 9: 0.0032649254 0.8412843 Trmt112,Wdr43,Riok2,Dimt1,Utp15,Riok1,...
## 10: 0.0041799164 0.8412843 Slc17a7,Adora1
## 11: 0.0042602556 0.8412843 Slc17a7,Adora1
## 12: 0.0046399072 0.8412843 Slc17a7,Ace2
## 13: 0.0049199016 0.8412843 Slc17a7,Ace2
## 14: 0.0053398932 0.8412843 Slc17a7,Ace2
## 15: 0.0054798904 0.8412843 Slc17a7,Ace2
## 16: 0.0057798844 0.8412843 Slc17a7,Folr1
## 17: 0.0062421973 0.8412843 Tcerg1,Gtf2f1,Brd4,Eapp,Cdk13,Cdk12,...
## 18: 0.0064102564 0.8412843 Trmt112,Wdr43,Kat2b,Riok2,Dimt1,Utp15,...
## 19: 0.0064401150 0.8412843 Cybb,Kcnmb1
## 20: 0.0077198456 0.8412843 Slc17a7,Folr1
## 21: 0.0091099498 0.8412843 Irf6,Erbb4
## 22: 0.0117797644 0.8412843 Cybb,Ncf2
## 23: 0.0121797564 0.8412843 Cybb,Ncf2
## 24: 0.0140800000 0.8412843 Cybb,Ncf2
## 25: 0.0152873379 0.8412843 Chrna7,Rims1
## 26: 0.0152996940 0.8412843 Slc17a7,Slc6a4,Grm4,Gpm6b
## 27: 0.0153159672 0.8412843 Slc6a4,Chrnb2
## 28: 0.0153367324 0.8412843 Agbl1,Irgm2
## 29: 0.0153917910 0.8412843 Entpd7,Ak3,Nme5,Nme4,Nme3,Nme1,...
## 30: 0.0154338843 0.8412843 Grm4,Grm8
## 31: 0.0154846454 0.8412843 Cybb,Ncf2
## 32: 0.0157155008 0.8412843 Cblc,Zgpat
## 33: 0.0161157025 0.8412843 Slc6a4,Slc44a4
## 34: 0.0167596648 0.8412843 Scin,Notch1
## 35: 0.0187011576 0.8412843 Nrg3,Tbx3
## 36: 0.0188800000 0.8412843 Chrna7,Rims1
## 37: 0.0194196116 0.8412843 Slc17a7,Prnp
## 38: 0.0200084119 0.8412843 Pde4c,Adcy4
## 39: 0.0201795964 0.8412843 Tnk1,Epha1,Pdgfra
## 40: 0.0207595848 0.8412843 Slc17a7,Grm4,Ntrk1
## pval padj leadingEdge\(~\)
\(~\)
Differntial expression between conditions for each gene
Exploratory data visualization for some Nephron markers.
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Six2", alpha = 1, border = TRUE) + ggtitle("Six2")
grid.text("(Fig38)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Cited1", alpha = 1, border = TRUE) + ggtitle("Cited1")
grid.text("(Fig39)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Crym", alpha = 1, border = TRUE) + ggtitle("Crym")
grid.text("(Fig40)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "C1qtnf12", alpha = 1, border = TRUE) + ggtitle("C1qtnf12/Fam132a")
grid.text("(Fig41)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Wnt4", alpha = 1, border = TRUE) + ggtitle("Wnt4")
grid.text("(Fig42)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "S100g", alpha = 1, border = TRUE) + ggtitle("S100g")
grid.text("(Fig43)", 0.15, 0.85)
plotSmoothers(GAM.1to18, assays(GAM.1to18)$counts, gene = "Umod", alpha = 1, border = TRUE) + ggtitle("Umod")
grid.text("(Fig44)", 0.15, 0.85)
\(~\)
\(~\)
To test differential expression between conditions, we use the “conditionTest” function implemented in tradeSeq. We performed the global (considering lineage1 and 2 together) and pairwise (each lineage separately) statistical tests and set the log fold change threshold at 1.2.
condRes <- conditionTest(GAM.1to18, global = T, pairwise = T , l2fc = log2(1.2))
# log2(1) = 0
# log2(1.2) = 0.2630344
# log2(1.4) = 0.4854268\(~\)
\(~\)
condRes$padj <- p.adjust(condRes$pvalue, "BH")
condRes$padj_lineage1 <- p.adjust(condRes$pvalue_lineage1, "BH")
condRes$padj_lineage2 <- p.adjust(condRes$pvalue_lineage2, "BH")
mean(condRes$padj <= 0.05, na.rm = TRUE)## [1] 0.0009775702# DEGs between wild-type and Wnt11 conditions
DEGs_general = subset(condRes , padj < 0.05)
DEGs_general## waldStat df pvalue waldStat_lineage1 df_lineage1
## Eif2s3x 36.28846 8 1.554845e-05 16.28014 4
## Xist 1319.19793 16 0.000000e+00 831.88335 8
## Pbdc1 93.78951 16 4.986012e-13 47.60188 8
## Hddc3 155.91543 10 0.000000e+00 82.94289 5
## Crebzf 89.11504 4 0.000000e+00 48.78038 2
## Gm35082 106.07208 6 0.000000e+00 51.02320 3
## Hbb-bt 181.70547 16 0.000000e+00 143.34267 8
## Hbb-bs 129.20132 16 0.000000e+00 47.78361 8
## Hbb-y 117.07229 16 0.000000e+00 16.15757 8
## Hba-x 65.32296 16 6.481468e-08 33.36763 8
## Hba-a1 83.99710 16 3.142919e-11 34.58898 8
## Hba-a2 83.24145 16 4.312095e-11 35.13281 8
## Col1a1 60.26111 16 4.729297e-07 30.96917 8
## Meg3 48.33821 16 4.202527e-05 24.81855 8
## Kdm5d 65.38066 6 3.606893e-12 37.19822 3
## Eif2s3y 154.58927 6 0.000000e+00 78.62146 3
## Uty 39.44500 4 5.637265e-08 26.00062 3
## Ddx3y 197.02849 4 0.000000e+00 103.88385 2
## pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Eif2s3x 2.665397e-03 20.00832 4 4.975146e-04
## Xist 0.000000e+00 487.31458 8 0.000000e+00
## Pbdc1 1.177078e-07 46.18763 8 2.189526e-07
## Hddc3 2.220446e-16 72.97254 5 2.464695e-14
## Crebzf 2.555489e-11 40.33466 2 1.743571e-09
## Gm35082 4.836731e-11 55.47784 4 2.579947e-11
## Hbb-bt 0.000000e+00 38.36279 8 6.450895e-06
## Hbb-bs 1.086656e-07 81.41771 8 2.531308e-14
## Hbb-y 4.017956e-02 100.91472 8 0.000000e+00
## Hba-x 5.286772e-05 31.95533 8 9.487318e-05
## Hba-a1 3.176323e-05 49.40812 8 5.309282e-08
## Hba-a2 2.528935e-05 48.10864 8 9.418119e-08
## Col1a1 1.422894e-04 29.29193 8 2.817783e-04
## Meg3 1.668499e-03 23.51966 8 2.757364e-03
## Kdm5d 4.177800e-08 28.18244 3 3.325465e-06
## Eif2s3y 1.110223e-16 75.96781 3 2.220446e-16
## Uty 9.534539e-06 13.79223 2 1.011707e-03
## Ddx3y 0.000000e+00 93.14465 2 0.000000e+00
## padj padj_lineage1 padj_lineage2
## Eif2s3x 1.684080e-02 1.000000e+00 5.376962e-01
## Xist 0.000000e+00 0.000000e+00 0.000000e+00
## Pbdc1 9.180743e-10 2.156878e-04 3.657106e-04
## Hddc3 0.000000e+00 8.137491e-13 7.751288e-11
## Crebzf 0.000000e+00 7.804464e-08 4.004330e-06
## Gm35082 0.000000e+00 1.266118e-07 6.771625e-08
## Hbb-bt 0.000000e+00 0.000000e+00 8.465878e-03
## Hbb-bs 0.000000e+00 2.156878e-04 7.751288e-11
## Hbb-y 0.000000e+00 1.000000e+00 0.000000e+00
## Hba-x 7.956218e-05 6.458320e-02 1.162070e-01
## Hba-a1 4.822548e-08 4.477150e-02 1.083860e-04
## Hba-a2 6.107585e-08 3.861683e-02 1.730391e-04
## Col1a1 5.442534e-04 1.533712e-01 3.235696e-01
## Meg3 4.298952e-02 1.000000e+00 1.000000e+00
## Kdm5d 6.037610e-09 9.569250e-05 4.699905e-03
## Eif2s3y 0.000000e+00 5.085932e-13 1.019906e-12
## Uty 7.414212e-05 1.588281e-02 9.783209e-01
## Ddx3y 0.000000e+00 0.000000e+00 0.000000e+00# DEGs between wild-type and Wnt11 conditions in lineage1
DEGs_lineage1 = subset(condRes , padj_lineage1 < 0.05)
DEGs_lineage1## waldStat df pvalue waldStat_lineage1 df_lineage1
## Xist 1319.19793 16 0.000000e+00 831.88335 8
## Pbdc1 93.78951 16 4.986012e-13 47.60188 8
## Hddc3 155.91543 10 0.000000e+00 82.94289 5
## Crebzf 89.11504 4 0.000000e+00 48.78038 2
## Gm35082 106.07208 6 0.000000e+00 51.02320 3
## Hbb-bt 181.70547 16 0.000000e+00 143.34267 8
## Hbb-bs 129.20132 16 0.000000e+00 47.78361 8
## Hba-a1 83.99710 16 3.142919e-11 34.58898 8
## Hba-a2 83.24145 16 4.312095e-11 35.13281 8
## Kdm5d 65.38066 6 3.606893e-12 37.19822 3
## Eif2s3y 154.58927 6 0.000000e+00 78.62146 3
## Uty 39.44500 4 5.637265e-08 26.00062 3
## Ddx3y 197.02849 4 0.000000e+00 103.88385 2
## pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist 0.000000e+00 487.31458 8 0.000000e+00
## Pbdc1 1.177078e-07 46.18763 8 2.189526e-07
## Hddc3 2.220446e-16 72.97254 5 2.464695e-14
## Crebzf 2.555489e-11 40.33466 2 1.743571e-09
## Gm35082 4.836731e-11 55.47784 4 2.579947e-11
## Hbb-bt 0.000000e+00 38.36279 8 6.450895e-06
## Hbb-bs 1.086656e-07 81.41771 8 2.531308e-14
## Hba-a1 3.176323e-05 49.40812 8 5.309282e-08
## Hba-a2 2.528935e-05 48.10864 8 9.418119e-08
## Kdm5d 4.177800e-08 28.18244 3 3.325465e-06
## Eif2s3y 1.110223e-16 75.96781 3 2.220446e-16
## Uty 9.534539e-06 13.79223 2 1.011707e-03
## Ddx3y 0.000000e+00 93.14465 2 0.000000e+00
## padj padj_lineage1 padj_lineage2
## Xist 0.000000e+00 0.000000e+00 0.000000e+00
## Pbdc1 9.180743e-10 2.156878e-04 3.657106e-04
## Hddc3 0.000000e+00 8.137491e-13 7.751288e-11
## Crebzf 0.000000e+00 7.804464e-08 4.004330e-06
## Gm35082 0.000000e+00 1.266118e-07 6.771625e-08
## Hbb-bt 0.000000e+00 0.000000e+00 8.465878e-03
## Hbb-bs 0.000000e+00 2.156878e-04 7.751288e-11
## Hba-a1 4.822548e-08 4.477150e-02 1.083860e-04
## Hba-a2 6.107585e-08 3.861683e-02 1.730391e-04
## Kdm5d 6.037610e-09 9.569250e-05 4.699905e-03
## Eif2s3y 0.000000e+00 5.085932e-13 1.019906e-12
## Uty 7.414212e-05 1.588281e-02 9.783209e-01
## Ddx3y 0.000000e+00 0.000000e+00 0.000000e+00# DEGs between wild-type and Wnt11 conditions in lineage2
DEGs_lineage2 = subset(condRes , padj_lineage2 < 0.05)
DEGs_lineage2## waldStat df pvalue waldStat_lineage1 df_lineage1
## Xist 1319.19793 16 0.000000e+00 831.88335345 8
## Pbdc1 93.78951 16 4.986012e-13 47.60188370 8
## Hddc3 155.91543 10 0.000000e+00 82.94288767 5
## Crebzf 89.11504 4 0.000000e+00 48.78038163 2
## Gm35082 106.07208 6 0.000000e+00 51.02319583 3
## Hbb-bt 181.70547 16 0.000000e+00 143.34267459 8
## Hbb-bs 129.20132 16 0.000000e+00 47.78361213 8
## Hbb-y 117.07229 16 0.000000e+00 16.15756801 8
## Nktr 28.70845 6 6.907214e-05 0.04676395 3
## Hba-a1 83.99710 16 3.142919e-11 34.58897703 8
## Hba-a2 83.24145 16 4.312095e-11 35.13280800 8
## Kdm5d 65.38066 6 3.606893e-12 37.19821884 3
## Eif2s3y 154.58927 6 0.000000e+00 78.62145884 3
## Ddx3y 197.02849 4 0.000000e+00 103.88384644 2
## pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist 0.000000e+00 487.31458 8 0.000000e+00
## Pbdc1 1.177078e-07 46.18763 8 2.189526e-07
## Hddc3 2.220446e-16 72.97254 5 2.464695e-14
## Crebzf 2.555489e-11 40.33466 2 1.743571e-09
## Gm35082 4.836731e-11 55.47784 4 2.579947e-11
## Hbb-bt 0.000000e+00 38.36279 8 6.450895e-06
## Hbb-bs 1.086656e-07 81.41771 8 2.531308e-14
## Hbb-y 4.017956e-02 100.91472 8 0.000000e+00
## Nktr 9.973478e-01 28.66169 3 2.637626e-06
## Hba-a1 3.176323e-05 49.40812 8 5.309282e-08
## Hba-a2 2.528935e-05 48.10864 8 9.418119e-08
## Kdm5d 4.177800e-08 28.18244 3 3.325465e-06
## Eif2s3y 1.110223e-16 75.96781 3 2.220446e-16
## Ddx3y 0.000000e+00 93.14465 2 0.000000e+00
## padj padj_lineage1 padj_lineage2
## Xist 0.000000e+00 0.000000e+00 0.000000e+00
## Pbdc1 9.180743e-10 2.156878e-04 3.657106e-04
## Hddc3 0.000000e+00 8.137491e-13 7.751288e-11
## Crebzf 0.000000e+00 7.804464e-08 4.004330e-06
## Gm35082 0.000000e+00 1.266118e-07 6.771625e-08
## Hbb-bt 0.000000e+00 0.000000e+00 8.465878e-03
## Hbb-bs 0.000000e+00 2.156878e-04 7.751288e-11
## Hbb-y 0.000000e+00 1.000000e+00 0.000000e+00
## Nktr 6.693817e-02 1.000000e+00 4.038426e-03
## Hba-a1 4.822548e-08 4.477150e-02 1.083860e-04
## Hba-a2 6.107585e-08 3.861683e-02 1.730391e-04
## Kdm5d 6.037610e-09 9.569250e-05 4.699905e-03
## Eif2s3y 0.000000e+00 5.085932e-13 1.019906e-12
## Ddx3y 0.000000e+00 0.000000e+00 0.000000e+00# Common DEGs between the two lineages
intersect(DEGs_lineage1 , DEGs_lineage2)## waldStat df pvalue waldStat_lineage1 df_lineage1
## Xist 1319.19793 16 0.000000e+00 831.88335 8
## Pbdc1 93.78951 16 4.986012e-13 47.60188 8
## Hddc3 155.91543 10 0.000000e+00 82.94289 5
## Crebzf 89.11504 4 0.000000e+00 48.78038 2
## Gm35082 106.07208 6 0.000000e+00 51.02320 3
## Hbb-bt 181.70547 16 0.000000e+00 143.34267 8
## Hbb-bs 129.20132 16 0.000000e+00 47.78361 8
## Hba-a1 83.99710 16 3.142919e-11 34.58898 8
## Hba-a2 83.24145 16 4.312095e-11 35.13281 8
## Kdm5d 65.38066 6 3.606893e-12 37.19822 3
## Eif2s3y 154.58927 6 0.000000e+00 78.62146 3
## Ddx3y 197.02849 4 0.000000e+00 103.88385 2
## pvalue_lineage1 waldStat_lineage2 df_lineage2 pvalue_lineage2
## Xist 0.000000e+00 487.31458 8 0.000000e+00
## Pbdc1 1.177078e-07 46.18763 8 2.189526e-07
## Hddc3 2.220446e-16 72.97254 5 2.464695e-14
## Crebzf 2.555489e-11 40.33466 2 1.743571e-09
## Gm35082 4.836731e-11 55.47784 4 2.579947e-11
## Hbb-bt 0.000000e+00 38.36279 8 6.450895e-06
## Hbb-bs 1.086656e-07 81.41771 8 2.531308e-14
## Hba-a1 3.176323e-05 49.40812 8 5.309282e-08
## Hba-a2 2.528935e-05 48.10864 8 9.418119e-08
## Kdm5d 4.177800e-08 28.18244 3 3.325465e-06
## Eif2s3y 1.110223e-16 75.96781 3 2.220446e-16
## Ddx3y 0.000000e+00 93.14465 2 0.000000e+00
## padj padj_lineage1 padj_lineage2
## Xist 0.000000e+00 0.000000e+00 0.000000e+00
## Pbdc1 9.180743e-10 2.156878e-04 3.657106e-04
## Hddc3 0.000000e+00 8.137491e-13 7.751288e-11
## Crebzf 0.000000e+00 7.804464e-08 4.004330e-06
## Gm35082 0.000000e+00 1.266118e-07 6.771625e-08
## Hbb-bt 0.000000e+00 0.000000e+00 8.465878e-03
## Hbb-bs 0.000000e+00 2.156878e-04 7.751288e-11
## Hba-a1 4.822548e-08 4.477150e-02 1.083860e-04
## Hba-a2 6.107585e-08 3.861683e-02 1.730391e-04
## Kdm5d 6.037610e-09 9.569250e-05 4.699905e-03
## Eif2s3y 0.000000e+00 5.085932e-13 1.019906e-12
## Ddx3y 0.000000e+00 0.000000e+00 0.000000e+00\(~\)
\(~\)
\(~\)
\(~\)
Bootstrap Hypothesis Test
First, we performed a t-test and Wilcoxon/Mann-whitney test on the distribution of pseudotime1 between wild-type and wnt11 groups. Both of which resulted in a p-value less than 0.05 percent.
# Comparing the means of two distribution with t test
t.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])##
## Welch Two Sample t-test
##
## data: slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## t = -5.1732, df = 1721.2, p-value = 2.57e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9289447 -0.4181948
## sample estimates:
## mean of x mean of y
## 3.274941 3.948511# Wilcoxon/Mann-whitney compares the medians of two distributions
wilcox.test(slingPseudotime(slingshot.wild.1to18)[, 1], slingPseudotime(slingshot.wnt11.1to18)[, 1])##
## Wilcoxon rank sum test with continuity correction
##
## data: slingPseudotime(slingshot.wild.1to18)[, 1] and slingPseudotime(slingshot.wnt11.1to18)[, 1]
## W = 322376, p-value = 1.517e-06
## alternative hypothesis: true location shift is not equal to 0\(~\)
\(~\)
Performing two Bootstrap hypothesis test based on mean and median of the distributions.
# In Bootstrap, we consider the null hypothesis is true, that is , the probability that
# a specific pseudotime comes from wild type is the same as it comes from Wnt11. Therefore,
# resampling is applied on aggregated data.
wild_pt1 = as.numeric(na.omit(slingPseudotime(slingshot.wild.1to18)[, 1]))
wnt11_pt1 = as.numeric(na.omit(slingPseudotime(slingshot.wnt11.1to18)[, 1]))
wild_pt1_mean = mean(wild_pt1)
wnt11_pt1_mean = mean(wnt11_pt1)
wild_pt1_median = median(wild_pt1)
wnt11_pt1_median = median(wnt11_pt1)
wild_pt1_mean## [1] 3.274941wnt11_pt1_mean## [1] 3.948511wild_pt1_median## [1] 2.841504wnt11_pt1_median## [1] 3.375997test_stat_mean = abs(wild_pt1_mean - wnt11_pt1_mean)
test_stat_median = abs(wild_pt1_median - wnt11_pt1_median)
L = length(c(wild_pt1, wnt11_pt1))
L## [1] 1732# Number of resampling
B = 10000
numbers = c(wild_pt1, wnt11_pt1)
# resampling L*B times
Bootstrap_matrix = matrix(sample(x = numbers , size = L*B , replace = T) , nrow = L, ncol = B)
dim(Bootstrap_matrix)## [1] 1732 10000Boot_stat_mean = rep(0,B)
Boot_stat_median = rep(0,B)
length(wild_pt1)## [1] 792for(i in 1:B){
Boot_stat_mean[i] = abs(mean(Bootstrap_matrix[1:792, i]) - mean(Bootstrap_matrix[793:1732, i]))
Boot_stat_median[i] = abs(median(Bootstrap_matrix[1:792, i]) - median(Bootstrap_matrix[793:1732, i]))
}
# Density plot of Bootstrap statistics
Boot_stat_mean[1:20]## [1] 0.08763163 0.08518836 0.15115668 0.03502281 0.05215037 0.14190328
## [7] 0.17647074 0.02751870 0.06701882 0.08203811 0.13314141 0.05520239
## [13] 0.01391266 0.20988616 0.22387382 0.02746928 0.13939191 0.05610408
## [19] 0.18180438 0.03027709dens = density(Boot_stat_mean)
plot(dens, frame = F , col = "steelblue", main = "Bootstrap statistics means")
polygon(dens, col = "steelblue")
grid.text("(Fig45)", 0.15, 0.75)
Boot_stat_median[1:20]## [1] 0.03267741 0.32980174 0.25902104 0.12543396 0.09649757 0.03764550
## [7] 0.50266934 0.06891447 0.11483682 0.25922714 0.24215826 0.12469570
## [13] 0.04605807 0.20465785 0.23964808 0.04441800 0.12099561 0.18582881
## [19] 0.24390804 0.09723747dens = density(Boot_stat_median)
plot(dens, frame = F , col = "indianred1", main = "Bootstrap statistics medians")
polygon(dens, col = "indianred1")
grid.text("(Fig46)", 0.15, 0.75)
# Computing the p-value for mean statistics and median statistics.
# Both Hypothesis tests gave rise to a p-value less than 0.05 percent.
sum(Boot_stat_mean > test_stat_mean)## [1] 0sum(Boot_stat_median > test_stat_median)/10000## [1] 0.0329mean(Boot_stat_median > test_stat_median)## [1] 0.0329\(~\)
\(~\)
Differential Expression Between Clusters
Differential expression between Wnt11KO vs Wild-type in all clusters. 1_wnt11-hom: wnt11 cells in cluster1 and so on.
so.small$Cluster_Genotype = paste(Idents(so.small), so.small$genotype, sep = "_")
Idents(so.small) = "Cluster_Genotype"
DEGs = FindMarkers(so.small, ident.1 = "1_wnt11-hom", ident.2 = "1_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 1.178851e-162 -3.3946361 0.056 1.000 2.186650e-158
## Gm35082 2.958657e-58 0.7788606 0.458 0.012 5.488013e-54
## Ddx3y 6.301919e-61 0.7426733 0.500 0.027 1.168943e-56
## Eif2s3y 9.122415e-53 0.6252948 0.432 0.017 1.692117e-48
## Ftl1-ps1 1.113321e-10 -0.5785791 0.596 0.762 2.065100e-06DEGs = FindMarkers(so.small, ident.1 = "21_wnt11-hom", ident.2 = "21_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 1.019964e-09 -3.990642 0 1 1.891931e-05DEGs = FindMarkers(so.small, ident.1 = "5_wnt11-hom", ident.2 = "5_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 2.186199e-100 -3.8232985 0.007 0.995 4.055181e-96
## Ddx3y 5.854359e-37 0.8172825 0.586 0.010 1.085925e-32
## Crabp1 2.507432e-08 0.8153917 0.500 0.285 4.651036e-04
## Cdkn1c 1.169589e-06 0.6940403 0.967 0.938 2.169471e-02
## Eif2s3y 3.342528e-29 0.6736892 0.516 0.026 6.200055e-25
## Gm35082 1.012469e-28 0.6589631 0.467 0.000 1.878028e-24
## Meg3 1.493468e-06 0.5397005 0.556 0.326 2.770234e-02
## Hddc3 2.344082e-15 0.5089163 0.539 0.192 4.348038e-11#subset(DEGs , abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)
DEGs = FindMarkers(so.small, ident.1 = "13_wnt11-hom", ident.2 = "13_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 1.586450e-38 -3.1064394 0.033 0.986 2.942706e-34
## Ftl1-ps1 1.476746e-13 -1.6014272 0.388 0.838 2.739216e-09
## Ddx3y 4.219953e-12 0.7150880 0.512 0.027 7.827591e-08
## Gm35082 7.515481e-11 0.6509924 0.430 0.000 1.394047e-06
## Eif2s3y 5.504918e-10 0.5853334 0.421 0.014 1.021107e-05DEGs = FindMarkers(so.small, ident.1 = "18_wnt11-hom", ident.2 = "18_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 9.585284e-18 -3.058509 0.043 1 1.777974e-13DEGs = FindMarkers(so.small, ident.1 = "15_wnt11-hom", ident.2 = "15_wild-type", verbose = F)
DEGs %>% arrange(desc(abs(avg_log2FC))) %>% filter(abs(avg_log2FC) > 0.5 & p_val_adj < 0.05)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## Xist 9.231529e-23 -2.971087 0.05 1 1.712356e-18